#### Set-Up ####
rm(list = ls())
Sys.info()[7]

# condition_data <- 0
# use_full_sample <- 0
expand <- 1
K <- 1000 # number of times to expand
N_bootstrap <- 1000
pa <- 0.005
pfp <- 0.01
pfn <- 0.01
J <- 1 

# Set file paths 
DATA <- "../external_links/real"
TEMP <- "../temp"
RESULTS <- "../output"
MODEL_RESULTS <- "../external/estimate_w_xs" # for separate dir

# Go to code folder
setwd("./")
sink("bootstrap_cf_compile.txt")
folder <- paste0(TEMP, "/cf_bs_iters")
if (file.exists(folder)) {
  cat("The folder already exists")
} else {
  dir.create(folder)
}
Sys.info()[7]
Sys.time()

# Packages
# library(MASS)
# library(tmvtnorm)
# library(haven)
# library(dplyr)
# library(tidyr)
# library(nloptr)
# library(corpcor)
library(stargazer)
# library(pracma)
# library(foreign)
# library(writexl)
# library(broom)
# library(png)
# library(tidyverse)

###################################################################################################################################################################
print("Begin time")
start_time <- Sys.time()
print(start_time)


for (b in 1:5) {
  bootstrap_iters_cf0_batch <- read.csv(paste0(TEMP, "/cf_bs_iters/bootstrap_iters_batch", b, "_cf0.csv"), header=FALSE)
  bootstrap_iters_cf1_batch <- read.csv(paste0(TEMP, "/cf_bs_iters/bootstrap_iters_batch", b, "_cf1.csv"), header=FALSE)
  bootstrap_iters_cf2_batch <- read.csv(paste0(TEMP, "/cf_bs_iters/bootstrap_iters_batch", b, "_cf2.csv"), header=FALSE)
  bootstrap_iters_cf4_batch <- read.csv(paste0(TEMP, "/cf_bs_iters/bootstrap_iters_batch", b, "_cf4.csv"), header=FALSE)
  bootstrap_iters_cf5_batch <- read.csv(paste0(TEMP, "/cf_bs_iters/bootstrap_iters_batch", b, "_cf5.csv"), header=FALSE)
  bootstrap_iters_cf7_batch <- read.csv(paste0(TEMP, "/cf_bs_iters/bootstrap_iters_batch", b, "_cf7.csv"), header=FALSE)
  bootstrap_iters_cf8_batch <- read.csv(paste0(TEMP, "/cf_bs_iters/bootstrap_iters_batch", b, "_cf8.csv"), header=FALSE)
  bootstrap_iters_cf9_batch <- read.csv(paste0(TEMP, "/cf_bs_iters/bootstrap_iters_batch", b, "_cf9.csv"), header=FALSE)

  if (b == 1) {
    bootstrap_iters_cf0 <- bootstrap_iters_cf0_batch
    bootstrap_iters_cf1 <- bootstrap_iters_cf1_batch
    bootstrap_iters_cf2 <- bootstrap_iters_cf2_batch
    bootstrap_iters_cf4 <- bootstrap_iters_cf4_batch
    bootstrap_iters_cf5 <- bootstrap_iters_cf5_batch
    bootstrap_iters_cf7 <- bootstrap_iters_cf7_batch
    bootstrap_iters_cf8 <- bootstrap_iters_cf8_batch
    bootstrap_iters_cf9 <- bootstrap_iters_cf9_batch
  } else {
    bootstrap_iters_cf0 <- rbind(bootstrap_iters_cf0, bootstrap_iters_cf0_batch)
    bootstrap_iters_cf1 <- rbind(bootstrap_iters_cf1, bootstrap_iters_cf1_batch)
    bootstrap_iters_cf2 <- rbind(bootstrap_iters_cf2, bootstrap_iters_cf2_batch)
    bootstrap_iters_cf4 <- rbind(bootstrap_iters_cf4, bootstrap_iters_cf4_batch)
    bootstrap_iters_cf5 <- rbind(bootstrap_iters_cf5, bootstrap_iters_cf5_batch)
    bootstrap_iters_cf7 <- rbind(bootstrap_iters_cf7, bootstrap_iters_cf7_batch)
    bootstrap_iters_cf8 <- rbind(bootstrap_iters_cf8, bootstrap_iters_cf8_batch)
    bootstrap_iters_cf9 <- rbind(bootstrap_iters_cf9, bootstrap_iters_cf9_batch)
  }
}
print(dim(bootstrap_iters_cf0))


sd0 <- apply(bootstrap_iters_cf0, 2, sd)
sd1 <- apply(bootstrap_iters_cf1, 2, sd)
sd2 <- apply(bootstrap_iters_cf2, 2, sd)
sd4 <- apply(bootstrap_iters_cf4, 2, sd)
sd5 <- apply(bootstrap_iters_cf5, 2, sd)
sd7 <- apply(bootstrap_iters_cf7, 2, sd)
sd8 <- apply(bootstrap_iters_cf8, 2, sd)
sd9 <- apply(bootstrap_iters_cf9, 2, sd)
out_sd <- rbind(sd0, sd1, sd2, sd4, sd5, sd7, sd8, sd9)
rownames(out_sd) <- c("Model", "First-best", "No testing", "Invasive only", "Free NIPT", "Out-of-pocket NIPT", "[1:51-1:200]", " > 1:201")
colnames(out_sd) <- c(c("Any testing", "NIPT only", "Invasive only", "NIPT followed by invasive", 
                        "Live birth", "live, chrom ab", "live, chrom ab, a_i > c_i", "live, chrom ab, a_i < c_i", "live, no chrom ab", 
                        "Terminated", "term, chrom ab, a_i > c_i", "term, no chrom ab", 
                        "Bad outcome", 
                        "G cost (p.c.)", "G cost, NIPT", "G cost, Inv", "Consumer surplus", "Share NIPT oop"))
print(paste0("number of bootstrap iterations: ", N_bootstrap))
print(t(out_sd))
stargazer(t(out_sd), digits = 2)
write.csv(out_sd, file=paste0(RESULTS, "/cf_outcomes_se.csv"), na=".", row.names=TRUE)

print("End time")
end_time <- Sys.time()
print(end_time)
time_diff <- end_time - start_time
print(time_diff)

#### EXIT ####

sink()
#rm(list = ls())